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ABSTRACT 

In  this  study,  an  optimization-based  approach  for 
simulating  the  walking  motion  of  a  digital  human  model  is 
presented.  A  spatial  skeletal  model  having  55  degrees  of 
freedom  is  used  to  demonstrate  the  approach.  Design 
variables  are  the  joint  angle  profiles.  Walking  motion  is 
generated  by  minimizing  the  mechanical  energy 
subjected  to  basic  physical  and  kinematical  constraints. 
A  formulation  for  symmetric  and  periodic  normal  walking 
is  developed  and  results  are  presented.  Backpack  and 
ground  reaction  forces  are  taken  into  account  in  the 
current  formulation,  and  the  effects  of  the  backpack  on 
normal  walking  are  discussed. 

INTRODUCTION 

Simulation  of  3D  human  walking  is  a  challenging 
problem  from  analytical  and  computational  points  of 
view.  An  accurate  biomechanical  walking  model  is 
needed  for  product  prototype  design  in  industry.  Several 
attempts  have  been  made  in  the  literature  to  develop 
realistic  mechanical  models.  One  approach  uses  motion 
capture  data  to  perform  simulations,  or  created  libraries 
that  can  be  used  to  query  for  a  particular  motion.  Motion 
capture  is  an  experiment-based  approach  where  the 
motion  of  a  subject  is  recorded  by  identifying  the  marker 
positions  in  the  Cartesian  coordinates.  The  joint  angle 
trajectories  are  then  generated  by  using  the  recorded 
data  and  the  inverse  kinematics.  This  approach  is  limited 
by  the  accuracy  of  the  experimental  data,  and  the 
simulated  motion  is  subject-specific  (McGuan,  2001). 

A  second  approach  is  the  ZMP-based  trajectory 
generation  method,  which  enforces  stability  of  the 
mechanism  during  walking.  With  the  preplanned  zero 
moment  point  (ZMP)  and  feet  location,  the  walking 
motion  can  be  generated  to  follow  the  desired  ZMP 
trajectory  using  an  optimal  control  strategy  (Yamaguchi 
et  al.,  1999;  Kajita  et  al.,  2003).  The  ZMP  concept  can 
also  be  incorporated  in  an  optimization  formulation  to 
synthesize  a  walking  pattern  by  maximizing  the  stability 
subject  to  physical  constraints  (Huang  et  al.,  2001 ;  Mu 


and  Wu,  2003;  Kim  et  al.,  2005).  The  inverted  pendulum 
model  is  another  method  that  has  been  used  to  solve  the 
walking  problem  because  biped  walking  can  be  treated 
as  an  inverted  pendulum.  Advantages  of  this  method  are 
its  simplicity  and  fast  solvable  dynamics  equations  (Park 
and  Kim,  1998).  However,  it  also  suffers  from  an 
inadequate  dynamics  model  that  cannot  generate  natural 
and  realistic  human  motion. 

In  addition  to  the  previously  mentioned  methods, 
optimization-based  trajectory  generation  is  aimed  at 
more  realistic  and  natural  humanoid  motion.  A  model 
with  large  degrees  of  freedom  (DOF)  is  not  a  problem, 
and  many  human-featured  criteria  can  be  considered 
simultaneously.  For  digital  human  simulations,  the 
objective  functions  represent  human  performance 
measures,  and  optimization  methods  are  used  to  solve 
for  the  feasible  joint  motion  profiles  that  optimize  an 
objective  function  and  satisfy  the  necessary  constraints 
(Chung  et  al.,  2007;  Thelen  et  al.,  2003).  Chevallereau 
and  Aousin  (2001)  planned  a  robotic  walking  and  running 
motion  using  optimization  to  determine  the  coefficients  of 
a  polynomial  approximation  for  profiles  of  the  pelvis 
translations  and  joint  angle  rotations.  Walking  was 
treated  as  a  combination  of  successive  single  support 
phases  with  instantaneous  double  support  phases 
defined  by  passive  impact.  Saidouni  and  Bessonnet 
(2003)  solved  for  cyclic,  symmetric  gait  motion  of  a  nine- 
DOF  model  that  moves  in  the  sagittal  plane.  The  control 
points  for  the  B-spline  curves  along  with  the  time 
durations  for  the  gait  stages  are  optimized  to  minimize 
the  actuating  torque  energy.  By  adopting  the  time 
durations  as  design  variables,  the  motion  for  the  single 
support  and  the  motion  for  the  double  support  are 
simultaneously  optimized.  Anderson  and  Pandy  (2001) 
developed  a  musculoskeletal  model  with  23  DOF  and  54 
muscles  for  normal  symmetric  walking  on  level  ground 
using  dynamic  optimization  method.  Muscle  forces  were 
treated  as  design  variables,  and  metabolic  energy 
expenditure  per  unit  distance  was  minimized.  The 
forward  dynamics  was  solved  for  kinematics,  and  ground 
reaction  forces  were  obtained  from  the  equilibrium 
condition  in  each  iteration.  The  repeatable  initial  and  final 
postures  were  given  from  the  experimental  data. 
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In  the  present  work,  an  optimization-based  predictive 
formulation  based  on  the  physics  of  motion  (the 
dynamics  of  the  motion)  is  developed;  it  relies  on  a  few 
experimental  data  or  not  at  all.  Instead  an  objective 
function  and  appropriate  constraints  are  defined  to 
simulate  natural  human  walking  motion.  This  is  called  the 
predictive  dynamics  approach.  The  digital  human  model 
consists  of  55  DOF,  where  a  DOF  characterizes  a  jointed 
pair  of  links  in  the  kinematics  sense,  where  various  links 
of  the  body  are  assumed  to  be  connected  by  revolute 
joints.  Six  DOF  represent  global  translation  and  rotation 
and  the  other  49  represent  the  kinematics  of  the  body. 
The  resultant  action  of  all  the  muscles  at  a  joint  is 
represented  by  the  torque  for  each  degree  of  freedom. 
The  torques  and  angles  at  a  joint  are  treated  as 
unknowns  in  the  optimization  problem.  The  cubic  13- 
spline  interpolation  is  used  for  time  discretization,  and 
the  well-established  robotic  formulation  of  the  Denavit- 
Flartenberg  method  is  used  for  kinematic  analysis  of  the 
mechanical  system.  The  recursive  Lagrangian 
formulation  is  used  to  develop  the  equations  of  motion 
and  was  chosen  because  of  its  known  computational 
efficiency.  The  approach  is  also  suitable  for  evaluation 
of  the  gradients  in  closed  form  that  are  needed  in  the 
optimization  process.  The  problem  is  formulated  as  a 
nonlinear  optimization  problem.  A  unique  feature  of  the 
formulation,  perhaps  the  most  significant  of  this  method, 
is  that  the  equations  of  motion  are  not  integrated 
explicitly.  They  are  imposed  simply  as  equality 
constraints  in  the  optimization  process,  thus  enforcing 
the  laws  of  physics.  A  program  based  on  a  sequential 
quadratic  programming  approach  is  used  to  solve  the 
nonlinear  optimization  problem.  The  control  points  for  the 
joint  angle  profiles  are  treated  as  design  variables.  For 
the  performance  measure,  the  mechanical  energy  that  is 
represented  as  the  integral  of  the  sum  of  the  squares  of 
all  the  joint  torques  is  minimized.  The  dynamic  stability  is 
achieved  by  satisfying  the  ZMP  constraint  throughout  the 
walking  motion.  The  solution  is  simulated  in  the  Santos™ 
environment  and  the  results  are  validated  and  verified  by 
examining  two  of  the  six  standard  determinants  from 
experiments  that  characterize  normal  walking.  In  addition 
to  the  normal  walking  case,  other  cases  that  simulate 
walking  with  a  shoulder  backpack  of  varying  loads  are 
also  addressed  and  the  results  are  demonstrated. 

SPATIAL  HUMAN  SKELETAL  MODEL 

The  kinematics  of  the  spatial  human  skeletal  model  of 
the  current  work  is  based  on  the  Denavit-Flartenberg 
method. 

DENAVIT-HARTENBERG  METHOD 

The  kinematics  relation  of  a  spatial  human  model  is 
represented  by  the  Denavit-Hartenberg  (DH)  method,  in 
which  4x4  homogeneous  transformation  matrices  relate 
two  adjacent  coordinate  systems  (Denavit  and 
Hartenberg,  1955).  The  DH  transformation  matrix 
includes  rotation  and  translation  and  is  a  function  of  four 
parameters:  6*  ,  dt ,  at ,  and  at ,  as  shown  in  Equation 
(1). 
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In  order  to  obtain  a  systematic  representation  of  a  serial 
kinematics  chain,  qe  Rn  is  defined  as  the  vector  of  n- 
generalized  coordinates,  the  joint  angles.  The  position 
vector  of  a  point  of  interest  in  the  Cartesian  space  can  be 
written  in  terms  of  the  joint  variables  as  X  =  X(q) .  In  this 

form,  the  augmented  4x1  vectors  °r„  and  rn  are 
defined  using  the  global  Cartesian  vector  X(q)  and  the 
local  Cartesian  vector  Xn  as: 


X(q)" 

;  = 

X' 

1 

5  n 
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where  X„  is  the  position  of  a  point  with  respect  to  the  nth 
coordinate  system.  Using  these  vectors,  °r„  can  be 
related  to  r„  as  follows: 

°rn  =  X  (q)  r„  (3) 

where 

°T„  (q)  =  H  =  °T1  (ch  )%{q2)-  n~lTn  (qn )  (4) 

i= 1 

55-DOF  WHOLE-BODY  MODEL 

A  spatial  digital  human  skeletal  model  with  55  DOF,  as 
shown  in  Figure  1 ,  is  considered  in  this  work.  The  model 
consists  of  six  physical  branches  and  one  virtual  branch. 
The  physical  branches  include  the  right  leg,  the  left  leg, 
the  spine,  the  right  arm,  the  left  arm,  and  the  head.  In 
these  branches,  the  right  leg,  the  left  leg,  and  the  spine 
start  from  the  pelvis,  while  the  right  arm,  left  arm,  and 
head  start  from  the  ending  joint  of  the  spine  ( z30 ,  z31 , 
z32  )■  The  virtual  branch  contains  six  global  DOFs, 
including  three  global  translations  ( z, ,  z2,  z3)  and  three 
global  rotations  ( z4 ,  z5,  z6)  located  at  the  pelvis,  and 
move  the  model  from  the  origin  (o-xyz)  to  the  current 
pelvic  position  (z4,  z5,  z6). 


Figure  1 .  The  55-DOF  digital  human  model 

DYNAMICS  MODEL 

FORWARD  RECURSIVE  KINEMATICS 


=  aa  ; 


r,  =  B  r  ; 


=  C„r„ 


(8) 


where  rn  represents  the  augmented  local  coordinates  of 
the  point  in  the  n,h  coordinate  system. 

BACKWARD  RECURSIVE  DYNAMICS 

Based  on  forward  recursive  kinematics,  the  backward 
recursion  for  the  dynamic  analysis  is  accomplished  by 
defining  a  4x4  transformation  matrix  D;  and  4x1 
transformation  matrices  E,.,  F. ,  and  G, ,  as  follows. 


Given  the  mass  and  inertia  properties  of  each  link,  the 
external  force  ftT  =  [/v  fy  f,  o] ,  and  the  moment 

htT  =  hy  h_  o]  for  link  /c(l<£<«)  defined  in  the 
global  coordinate  system,  the  joint  actuation  torques  rf 
for  /  =  nto  1  are  computed  as  follows  (Toogood,  1989): 
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In  this  process,  4x4  transformation  matrices  A;  ,  By , 
and  Cj  are  defined  to  represent  the  recursive  position, 

velocity,  and  acceleration  for  the  /h  joint  respectively. 
Given  the  link  transformation  matrix  (  T.  )  and  the 

kinematics  state  variables  for  each  joint  ( q} ,  qt ,  and 
c 'jj ),  then  for  j  =1  to  n  we  have: 

A  =  T  T:  Ts  •  •  •  T,  =  A  ,  T  (5) 


G,  =ht8a  +Gi+1  (13) 

where  Dn+1  =  En+1  =  F„+1  =  Gn+1  =  [  0  ]  ;  J,  is  the  inertia 
matrix  for  link  /';  mi  is  the  mass  of  link  /';  g  is  the  gravity 
vector;  'r.  is  the  location  of  center  of  mass  of  link  /  in  the 
local  frame  /;  krf  is  the  position  of  the  external  force  in 

the  local  frame  k;  z0  =  [0  0  1  0]T  ,  and  8a  is 
Kronecker  delta. 
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A0  =  [  I  ]  and  B0  =  C0  =  [  0  ] . 

After  obtaining  all  transformation  matrices  A; ,  ,  and 

C j ,  the  global  position,  velocity,  and  acceleration  of  a 

point  in  the  Cartesian  coordinate  system  can  be 
calculated  using  the  following  formulas: 


The  first  term  in  the  torque  expression  is  the  inertia  and 
Coriolis  torque,  the  second  term  is  the  torque  due  to 
gravity,  the  third  term  is  the  torque  due  to  external  force, 
and  the  fourth  term  represents  the  torque  due  to  the 
external  moment. 

SENSITIVITY  ANALYSIS 

The  gradients  of  the  torque  for  the  3D  human 
mechanical  system  with  respect  to  the  state  variables 

—L- ,  — r1 - ,  -—7-  (/  =  1  to  n;  k  =  1  to  n)  can  be  evaluated  in 
oqk  oqk  oqk 

a  recursive  way  using  the  foregoing  recursive  Lagrangian 
dynamics  formulation  (Equations  (5)-(13)). 
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SYMMETRIC  NORMAL  GAIT 

In  the  current  work,  normal  walking  is  assumed  to  be 
symmetric  and  cyclic.  As  defined  in  the  literature,  a 
complete  gait  cycle  includes  two  continuous  steps  (one 
stride).  The  step  is  further  divided  into  two  phases,  the 
single  support  phase  and  the  double  support  phase. 

Considering  the  foot  flexion,  the  two  support  phases  can 
be  detailed  into  four  basic  supporting  modes:  right  foot 
single  support  (RSS),  left  foot  single  support  (LSS),  right 
foot  leading  double  support  (RDS)  and  left  foot  leading 
double  support  (LDS),  as  shown  in  Figure  2.  In  this  work, 
the  gait  cycle  starts  from  the  left  heel  strike  (LDS)  and 
goes  through  the  right  swing  (LSS),  the  right  heel  strike 
(RDS),  and  the  right  swing  (RSS)  before  coming  back  to 
left  heel  strike  (LDS)  again.  By  applying  the  symmetry 
conditions,  it  is  possible  to  consider  only  one  step  in  a 
normal  steady  gait  cycle  instead  of  considering  two 
steps.  As  a  result,  the  computational  costs  are 
significantly  minimized,  as  depicted  in  Figure  3. 
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Figure  2.  Four  basic  supporting  modes  (top  view) 


Figure  3.  A  normal  step  with  initial  and  final  postures 
symmetry 

GROUND  REACTION  FORCES 

In  the  single  support  phase,  one  foot  supports  the  whole 
body  and  the  ZMP  stays  in  the  foot  area  so  that  ground 
reaction  forces  (GRF)  can  be  applied  at  the  ZMP  directly. 
However,  in  the  double  support  phase,  the  ZMP  is 
located  between  the  two  supporting  feet,  and  therefore, 
the  resultant  GRF  needs  to  be  distributed  into  the  two 
feet  appropriately.  This  distribution  process  can  be 
treated  as  a  sub-optimization  problem  (Dasgupta  and 
Nakamura,  1999).  In  order  to  simplify  this  process,  a 
linear  relationship  is  used  to  distribute  GRF  in  the  double 
support  phase  as  shown  in  Figure  4  (a). 


My 


Figure  4.  Applying  ground  reaction  forces  on  ZMP 

(a)  double  support  phase,  (b)  single  support  phase 

In  Figure  4,  points  A  and  B  (triangles)  are  center  points 
of  the  middle  joints  of  the  feet.  dl  and  d2  are  the 

distances  from  the  ZMP  (circles)  to  points  A  and  B, 
respectively.  Note  that  there  are  only  normal  moment 
My  and  resultant  force  R  ( Rx ,  Ry ,  Rz )  at  ZMP,  as  M_ 

and  Mx  vanish  due  to  the  assumption  of  no-slip 

condition.  The  GRF  is  linearly  decomposed  to  the  central 
points  as  follows: 
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In  this  work,  a  two-step  algorithm  is  used  to  calculate 
GRF  as  depicted  in  the  following  flowchart. 


where  ndof  is  the  number  of  degrees  of  freedom. 
CONSTRAINTS 

Several  constraints  are  proposed  and  implemented  in 
this  work  to  satisfy  laws  of  physics  and  boundary 
conditions  through  out  the  normal  walking  process. 
These  constraints  include  Joint  angle  limits,  ground 
penetration,  foot  contacting  positions,  ZMP  stability, 
Pelvic  velocity,  soft  impact,  Knee  flexion  at  mid-swing, 
symmetry  conditions,  and  the  arm-leg  coupling. 

Joint  Angle  Limits 

The  joint  angle  limits  accounting  for  the  physical  range  of 
motion  are  obtained  from  experiments: 


Figure  5.  Flowchart  of  two-step  algorithm  to  calculate 
GRF  (a)  without  GRF,  and  (b)  with  GRF 

First,  given  current  state  variables,  the  recursive 
Lagrangian  dynamics  of  the  whole  body  is  used  to 
calculate  torques  without  GRF.  The  resulting  global 
forces  in  the  virtual  branch  are  not  zero  because  of 
excluding  GRF;  however,  the  moments  about  the  x  and  z 
at  ZMP  are  zero  due  to  the  ZMP  calculation. 

Second,  GRF  are  assumed  to  be  acting  at  the  ZMP. 
Considering  the  equilibrium  of  the  global  forces  and 
moments  in  the  virtual  branch  with  the  ground  reaction 
forces  and  moments,  three  forces  and  moment  about  the 
y  axis  at  the  ZMP  are  calculated.  The  GRF  are  then 
treated  as  external  forces  and  moments  for  the  entire 
model  and  the  updated  joint  torques  are  recovered  from 
the  equations  of  motion. 

OPTIMIZATION  FORMULATION 

DESIGN  VARIABLES 

In  the  current  formulation,  the  design  variables  are  the 
joint  profiles  qft)  for  a  symmetric  and  cyclic  gait  motion. 

Besides  the  joint  profiles,  the  initial  posture  is  also 
optimized  rather  than  specified  from  the  experiment. 
Also,  the  final  posture  should  satisfy  the  symmetry 
condition  with  the  initial  posture  so  that  continuous  joint 
profiles  are  generated.  Meanwhile,  the  torque  profiles  are 
calculated  from  joint  profiles  using  the  recursive 
Lagrangian  dynamics  equations. 

OBJECTIVE  FUNCTION 


qL<q(0<qc/,  0  <t<T  (19) 

where  qL  is  the  lower  limit  on  the  joint  angles  and  q77  is 
the  upper  limit  on  the  joint  angles. 

Ground  Penetration 


Bipedal  walking  is  characterized  with  unilateral  contact  in 
the  whole  process.  While  the  foot  is  in  contact  with  the 
ground,  the  height  of  the  contacting  points  is  zero.  The 
other  points  should  be  above  the  ground  and  the  height 
greater  than  zero  (Figure  6),  as  follows: 


yft)  =  0,  0  <t<T,  i e  contact 

yi (t)  >  0,  0 <t  <T,  it  contact 


(20) 


Figure  6.  Foot  ground  penetration  (circle:  on  the  ground; 
triangle:  above  the  ground) 


Contacting  Position 


The  energy-related  performance  measure,  the  integral  of 
the  squares  of  all  joint  torques,  is  used  as  the  objective 
function  for  the  walking  motion  prediction  that  is  defined 
as: 


Since  the  step  length  L  is  given,  the  foot  contacting 
position  is  known  and  specified  at  each  time. 

x;  (t)  =  x; ,  i  e  contact  (21) 


T  ndof 

Minimize  /(q)  =  |  f'  xf  (q )dt 

t= o  i=1 


(18) 


where  x,  is  the  specified  contacting  position. 
ZMP  Stability 


The  stability  is  achieved  by  constraining  the  ZMP  to  be  in 
the  foot  supporting  region  (FSR)  (Vukobratovic  and 
Borovac,  2004). 


The  arm-leg  coupling  motion  guarantees  that  the  swing 
of  the  arm  and  leg  on  the  same  side  are  in  the  opposite 
directions. 


zZMP(t)eFSR,  xZMP(t)e  FSR  0  <t<T  (22) 

Pelvic  Velocity 


bright  _  arm  bright Jegit)^  0» 

V-am,  (tHeft.leg  (0  ^  °,  0<t<T 


The  pelvic  translational  velocity  along  the  walking 
direction  in  normal  walking  is  stable  and  almost  a 
constant  ( V)  quantity;  this  is  treated  as  a  constraint. 

W')  =  V  0  <t<T  (23) 

Soft  Impact 

The  landing  impact  issue  is  taken  into  account  by 
imposing  velocity  of  contacting  points  to  be  zero.  This  is 
the  so-called  soft  impact,  which  makes  the  gait  motion 
gentle  and  smooth. 

x,(t)  =  0,  y,(f)  =  0,  z,(t)  =  0,  ^ 

0<t  <T,  i  e  contact 

Knee  Flexion  at  Mid-swing 

Knee  flexion  at  mid-swing  is  one  of  the  six  determinants 
of  normal  gait  (Ayyappa,  1997).  Biomechanical 
experiments  show  that  the  maximum  knee  flexion  of 
normal  gait  is  around  60  degrees  regardless  of  age  and 
gender.  In  addition,  this  also  ensures  enough  space 
between  the  foot  and  the  ground  to  avoid  collision.  This 
constraint  is  imposed  as 

Qknee  (0  =  60  +  8,  t  =  (25) 

where  e  is  a  small  range  of  motion. 

Symmetry  Condition 

The  gait  simulation  starts  from  the  left  heel  strike  and 
ends  with  the  right  heel  strike.  The  initial  and  final 
postures  should  satisfy  the  symmetry  conditions.  These 
conditions  are  expressed  as 


NUMERICAL  DISCRETIZATION 

For  the  optimization  problem,  the  entire  time  domain  is 
discretized  by  B-spline  curves,  which  are  defined  by  a  set 
of  control  points  P  and  time  grid  points  (knots)  t.  A  B- 
spline  is  a  numerical  interpolation  method  that  has  many 
important  properties,  such  as  continuity,  differentiability, 
and  local  control.  These  properties,  especially 
differentiability  and  local  control,  make  B-splines 
competent  to  represent  joint  angle  trajectories,  which 
require  smoothness  and  flexibility. 

Let  t  =  {  t0,  f[ , tm)  be  a  non-decreasing  sequence  of 

real  numbers,  i.e.,  t,  <ti+l,  i  =  0,..  ,m- 1.  The  tt  are  called 

knots,  and  they  are  non-decreasingly  spaced  for  a  B- 
spline.  A  cubic  B-spline  is  defined  as 

net 

9(0  =  2^(0^;  0<*<r  (28) 

j= o 

where  the  {/>},  j  =  0,..,nct  are  the  (nct+ 1)  control 

points,  and  the  {Nj4  (/)}  are  the  cubic  B-spline  basis 
functions  defined  on  the  non-decreasing  knot  vector. 

Since  the  first-  and  second-order  derivatives  of  the  joint 
angles  are  needed  in  the  optimization  problem,  the 
derivatives  of  a  cubic  B-spine  curve  can  be  easily 
obtained  because  only  the  basis  functions  are  functions 
of  time.  Therefore,  the  original  continuous  variable 
optimization  problem  is  transformed  into  a  parameterized 
optimization  problem  by  using  Equations  (28)  and  the 
following  two  equations: 

net 

i{t)  =  llNjA{t)Pj-,  0<t<T  (29) 

j= o 


diji-ft  (6)  —  <li_right  (T) 

qjJ0)=qjJT), 
qjY(0)  =  -q]Y(T), 
q  jZ  (0)  =  ~qjz  CO, 


(26) 


where  x,  y,  z  are  the  global  axes,  and  the  subscript  / 
represents  the  DOF  of  the  leg,  arm  and  shoulder  joints 
and  j  represents  other  DOF. 


net 

</(')  =  Z'V  ;(,)/>:  0 <t<T  (30) 

7=0 

q ,  q ,  and  q  are  functions  of  t  and  P;  therefore,  torque 
r  =  r(t,P)  is  an  explicit  function  of  the  knot  vector  and 
control  points  from  the  equation  of  motion.  Thus,  the 
derivatives  of  a  torque  t  with  respect  to  the  control  points 
is  computed  using  the  chain  rule  as 


Arm-Leg  Coupling 


dr 

dPt 


dz  dq  dz  dq  dz  dq 


dq  ()!)  dq  dPl  dq  dPl 


(31) 


Finally,  five  control  points  are  used  for  each  DOF  and 
thus  we  have  5x55  =  275  design  variables  and  629 
nonlinear  constraints. 

NUMERICAL  RESULTS 

A  large-scale  sequential  quadratic  programming  (SQP) 
approach  in  SNOPT  (Gill  et  al.,  2002)  is  used  to  solve 
the  nonlinear  optimization  problem  of  normal  walking.  To 
use  the  algorithm,  cost  and  constraint  functions  and  their 
gradients  need  to  be  calculated.  The  foregoing 
developed  recursive  kinematics  and  dynamics 
formulations  provide  accurate  gradients  to  improve  the 
computational  efficiency  of  the  numerical  optimization 
algorithm.  The  appropriate  normal  walking  parameters 
(velocity  and  step  length)  are  obtained  from 
biomechanics  literature  (Inman  et  al.,  1981).  In  addition 
to  normal  walking,  the  current  work  also  considered 
situations  where  people  walk  and  carry  backpacks  with 
various  weights  (0  lbs,  20  lbs,  and  40  lbs).  Each  case 
requires  about  600  cpu-seconds  on  a  Pentium(R)  4  3.46 
GHz  computer. 

NORMAL  GAIT  MOTION 

To  solve  a  normal  gait  motion,  the  user  inputs  a  pair  of 
parameters,  for  example,  normal  walking  velocity  V=  1.2 
m/s  and  step  length  L  =  0.7  m.  Figure  7  shows  the 
resulting  stick  diagram  of  a  3D  human  walking  on  flat 
ground.  As  expected,  correct  bending  of  knee  occurs  to 
avoid  collision,  and  the  arms  swing  to  balance  the  leg 
swing  and  GRF.  The  continuity  condition  is  satisfied  to 
generate  a  smooth  walking  motion  while  the  initial 
posture  is  also  optimized.  It  is  important  to  note  that  the 
spine  remains  upright  automatically  to  reduce  energy 
expenditure  in  the  walking  motion.  In  addition,  the  ZMP 
trajectory  is  smooth  and  stays  in  the  support  region 


Figure  7.  Stick  diagram  of  optimized  cyclic  walking 
motion  (square  dot  is  ZMP) 

Figure  8  depicts  the  torque  profiles  of  the  right  hip,  right 
knee,  and  right  ankle.  It  is  evident  that  the  ankle  torque  is 
quite  small  in  the  swing  phase  and  that  the  peak  value 
occurs  at  the  heel  strike. 


Figure  8.  Joint  torque  profiles  of  right  leg 

NORMAL  WALKING  WITH  BACKPACK 

In  this  process,  a  shoulder  backpack  is  considered  in  the 
walking  motion  prediction  with  the  walking  velocity  V  = 
1.0  m/s  and  step  length  L  =  0.5  m.  Three  cases  are  tried 
with  varying  backpack  weights:  0  lbs,  20  lbs,  and  40  lbs. 
3D  diagrams  of  the  walking  motion  are  compared  in 
Figure  9,  and  reasonable  spine  bending  is  observed.  The 
GRF  are  depicted  as  in  Figure  10,  and  a  heavier 
backpack  results  in  larger  GRF.  The  torque  profile  of  the 
right  knee  is  illustrated  in  Figure  1 1 . 


Figure  9.  Normal  walking  with  backpack:  (a)  0  lbs;  (b)  20 
lbs;  (c)  40  lbs 


Percent  of  Cycle 


Figure  10.  Vertical  GRF  of  walking  with  backpack 


Time  (s) 


Figure  1 1 .  Torque  profile  of  right  knee  with  backpack 

VALIDATION 

In  this  work,  two  predicted  walking  determinants,  the 
knee  and  ankle  joint  angle  profiles,  were  chosen  for 
comparison  purposes  with  data  obtained  from  four 
normal  subjects  as  shown  in  Figure  12.  These  walking 
parameters  are  a  subset  of  the  six  well-known  walking 
determinants  defined  by  Saunders  et  al.  (1953). 

The  data  predicted  by  the  dynamic  model  for  the  knee 
and  ankle  joints  were  also  added  to  Figure  12.  The 
graphs  clearly  show  that  the  predicted  trend  is  in  the 
range  of  normal  in  terms  of  joint  angle  amplitude  and 
time  history. 


Heel  Strike  0% 


Heel  Strike  1 00% 


Percent  of  Cycle 


Figure  12.  Verification  of  joint  angle  profiles  of  (a)  knee 
and  (b)  ankle  with  experiments 

CONCLUSION 

In  this  paper,  the  motion  prediction  of  the  normal  gait  of  a 
55-DOF  digital  human  model  was  presented.  Based  on 
the  experimental  data  and  biomechanical  requirements, 
the  results  have  demonstrated  the  ability  of  the  proposed 
methodology  to  predict  realistic  human  motions  with  the 
complex  spatial  skeleton  model.  The  normal  gait  was 
treated  as  a  cyclic  and  symmetric  motion  with  repeatable 
initial  and  final  postures.  The  motion  planning  was 
formulated  as  a  large-scale  nonlinear  programming 
problem.  Joint  profiles  were  discretized  using  cubic  B- 
splines,  and  the  corresponding  control  points  were 
treated  as  unknowns.  The  energy-related  objective 
function,  which  represents  the  integral  of  squares  of  all 
joint  torques,  was  minimized.  The  arm  motion  was 
incorporated  by  considering  the  role  of  the  yawing 
moment  and  the  arm-leg  coupling  motion  constraints  in 
the  formulation.  The  walking  determinants  were  obtained 
to  verify  the  gait  motion.  The  effect  of  backpack  on  gait 
motion  was  studied,  and  plausible  posture  responses 
were  achieved.  The  kinetic  data  such  as  joint  torque  and 
ground  reaction  forces  were  also  analyzed.  The  limitation 
of  current  gait  formulation  is  that  only  energy-saved 


motion  is  generated  and  the  walking  motion  is  on  a  level 
ground  with  symmetry  conditions.  Walking  on  uneven 
terrain  and  consideration  of  other  performance  measures 
such  as  maximizing  speed,  minimizing  fatigue  etc.  are 
ongoing  research  topics. 
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